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Abstract 

We consider a gravitational field in steady state galaxy models of two kinds. Some of 
them are axisymmetrical and others are triaxial. Equipotentials and potential law are given 
separately in accordance to Kutuzov and Ossipkov (1980a). The relatively simple potential 
law is based on Kuzmin-Malasidze model (1969). Two kinds of models contain four and five 
structural parameters respectively. One composite model is suggested as well. Some examples of 
trajectories are calculated in these models. The simplest method to describe orbits is drawing 
their projections on coordinate planes. However it needs a great amount of calculation and 
makes troubles in an interpretation of information. In the case of axisymmetrical models a 
motion in co-moving meridional plane (with cylindrical coordinates R, z) is considered as a 
common way. In the case of triaxial models one can use three different co-moving planes passing 
through moving star and corresponding coordinate axis. We describe models in sections 2-4, 
calculated orbits are discussed in section 5. 
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1 Introduction 

In recent years trajectories of stars in various models of the gravitational field of the galaxies have 
been constructed intensively. We intend to clear up how do orbit characteristics depend on model 
properties. Similar investigation was fulfilled by Kutuzov & Ossipkov (1992), where box eccentric- 
ities and maximal elevation of 70 orbits of open clusters in three models were compared. Here we 
investigate relation of the orbit properties with variation of the properties of different in principle 
models — axially symmetric and triaxial. This work continues our research (Kutuzov & Raspopova 
2008). 

There is a great variety of models in literature. We mention just few of them. Ossipkov (1997) 
has performed an exhaustive analysis of fourth-order (quaternary) equipotentials of general form 
with some parameters for systems with rotational symmetry. 

Dynamical models of triaxial elliptical galaxies were constructed by Arnold et al. (1994) using 
Stackel potential for calculating the dependence velocity field of the model on the shapes of orbits. 
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Merritt and Fridman (1996) showed essential role of the stochastic orbits in triaxial gravitational 
field with central density cusps. For similar model Merritt and Valluri (1996) found that central 
point mass designed to represent nuclear black hole causes stochastic orbits to diffuse through phase 
space. Kutuzov (1997) suggested the model with eighth-order equipotentials, providing triaxial mass 
distribution. 

Here we use simpler triaxial fourth-order equipotentials of special form (Kutuzov 1998). One of 
the parameters, responsible for the triaxial shape, is specified in the form of a variable in such a way 
that equipotentials asymptotically tend to spheres, thus, the natural requirement for the models of 
finite mass is satisfied. A rarely used property, according to which the dependence of the coefficients 
in the equation of level surfaces on the parameters of the family does not change the order of the 
equation, is employed. In the special case a biaxial disk embedded in the halo occurs. 

2 Family of axially symmetric models 

All the quantities are dimensionless. The change to dimensional quantities is fulfilled by multipli- 
cation ones by the corresponding scaling parameters. Basic ones are units of length, potential and 

mass 

f, $, M = r$/G, 

where G is gravitational constant. 

At first, we consider axisymmetrical models. The potential law is specified analytically with free 
parameters and the one independent variable £. The only argument £ is a function of coordinates 

¥> = </?(£) , e 2 = w, z) . (i) 

The variable £ is the equipotential parameter that changes from one equipotential to another, deter- 
mining their family. Putting £ = const formula (pQ) gives equation of a fixed equipotential. 

Equipotentials and potential law are given separately in accordance to Kutuzov & Ossipkov 
(1980a). Ten models being the special cases of mentioned family were listed by Kutuzov & Ossipkov 
(1980b). These models were suggested earlier by P. P. Parenago (1950, 1952), G. M. Idlis (1961), 
G. G. Kuzmin (1953, 1956), A. Toomre (1963), H. Plummer (1911), M. Miyamoto and R. Nagai 
(1975), M. Henon (1959), G. G. Kuzmin and U.-I. K. Veltmann (1972). Euipotentials, constructed 
there, were generalized later (Kutuzov 1989). According to ([I]), equipotentials are determined by the 
function f(R, z). Two-parametric family of equipotentials was suggested (Kutuzov 1989)as follows 

i 2 = r 2 + 2/2(1 -e)(q-e) , r 2 = R 2 + z 2 , q = Ve 2 + z 2 . (2) 

Quantities e, fj, are structural parameters of the model 

< e < 1 , 0</i<oo. 
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The first term on the right-hand of equation describes the spherical shape of equipotentials at 
infinity. The second term causes to a disk component in the model when e — > 0. These equipotentials 
coincide with the ones of the model of Miyamoto & Nagai (1975) for fi = 1. Let us compare our 
equipotentials with the model of Satoh (1980). We perform the dimensional potential of that model 
to 

HO = : J^r =a> e = r 2 + 2(l-e)(q-e). 
- e 2 + £ 2 



It is obvious that these equipotentials are a special case of our ones (|2j) when \i = 1. For e = 1 the 
model is spherical, potential has discontinuity in the center, where £ = 0. System degenerates into 
mass point with Kepler potential. 
We take another potential law 

¥>(£) = " ^==3 £ = (3) 

that coincides in the plane z = with Kuzmin-Malasidze law(1969). Thus accepted family of the 
models has four structural dimensionless parameters e, /i, a, x. It gives great opportunities for 
modeling of the various galaxies and star clusters. The potential is normed so that in the center 
<p(Q) = 1. 

Asymptotics of the potential allows us to find the expression for the dimensionless mass of the 
model 

M = 



When e = 1 or fi = the force field is spherically symmetrical, because £ = r according to In 
the case of e = the force has a discontinuity in the z = plane, as q = \z\. This determines the 
existence of an infinitely thin circular disk, imbedded into continuous halo; dimensionless mass of 
the disk is equal fi (Kutuzov 1989). 



3 Family of triaxial models 

We still take potential in the form of but now its argument £ is the function of three Cartesian 
coordinates x, y, z (with the origin at center of mass) 

£ 2 = g(x,y,z) . 

Clearly, the addition of £ to the list of arguments, i.e. the representation f(x, y,z,£), does not change 
the order of the equipotential equation relative to the x, y, z coordinates. 
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Let the function f(x,y,z,£) consist of three terms 



e = r 2 + 2 - £ )(q-e) + r(£) s 2 , < £ < oo . (4) 

Here. 



r 



2 = x 2 + y 2 + z 2 , q = Ve 2 + z 2 



2 _ 2 2 

s = y — x 



> 0, \x\ < \y\ 
< 0, \x\ > \y\ 



The new expression differs from (T5]) by his third term, that is responsible for the triaxial shape. 
Coefficient r is assumed to be a function of £. Let us call it the triaxiality parameter. Sometimes it 
is convenient to express r in terms of some bounded function x{0 that is proportial to it 

r(0 = (1 - e) X(0 ■ (5) 

For e = 1 r 2 remains on the right side of eq. (jl]), which implies that the model is spherical. But 
now equality /i = does not mean transfer to the spherical model, because the triaxiality parameter 
might be not equal zero. Below we'll discuss function r(£) in details. We only note here that for 
r = the model ceases to be triaxial, including models with rotational symmetry (T5]) as a special 
case. Note also that the parameter £ is expressed in terms of the x, y, z coordinates implicitly. To 
calculate it from the coordinates, we must solve eq. (jl]) using eq. ([5]), for example, by successive 
approximations by assuming initially that r = 0. 

The equation of equipotential (jl]) has fourth order with respect to coordinates. It gives the 
second-order equation in equatorial plane 

e = (l-r)x 2 + (l + r)y 2 . (6) 

The section of the equipotential by the equatorial plane is an ellipse for |r| < 1. Equation ([6]) yields 
the following expressions for the semi-major and semi- minor axes of the ellipse 

If ^ r(£) < 1, the two axes are real, and b ^ £ ^ a . 

Kutuzov (1998) suggested the following constraints on the triaxiality parameter t(£). It is de- 
sirable (but not necessary) that, as one approaches the system's center, the ellipse (EJ) changes to a 
circumference that would shrink to a point 

r(0) = . (8) 
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For the potential to become spherically symmetric at infinity, the following limit must exist 

lim r(0 = . (9) 

In this case, r must tend to zero as some negative power of £. The third term on the right side of 
(J3j) will then have an order that is less than 2. Since the second term has only the first order in z, 
the first term, which always has order 2, dominates. 

In order that the semimajor axis a of the ellipse ([6]) be no smaller than its semiminor axis b, and 
that both axes remain bounded for a finite £, we assume that 

< r(0 < 1 . (10) 

The equipotentials of the family must not intersect (and touch) each other. This requires that 
the axes a and b be monotonically increasing functions of £ 

a'(£) > , b'(0 > . 

The prime denotes differentiation with respect to £ 2 . Taking into account (J7J), we obtain differential 
inequalities which the triaxiality parameter must satisfy 

£V(£) - r(0 + 1 > 0, £V(0 - r(0 - K 0. 

Here are three examples of the function x(0> that satisfy the necessary conditions 

with the following constraint on the positive parameters 

^<-— , 1 = 1,2,3. 

1 — e 

The first function does not satisfy constraint jS]), so that the model is also triaxial in the central 
region. The function r is defined in terms of x by formula (jSJ). The coefficients allow us to make 
the triaxiality parameter r arbitrarily small. 

We still use expression (j3J for the potential law. In the case when the coefficients A4 are equal 
zero, the model coincides with axially symmetric one entirely. 

A circular velocity has no sence in triaxial model. So we introduce quasi-curcular velocity, that 
just characterizes the gravitational field 

, ^{R,z) = y{t;)\ x=y , R=^Tf. (12) 

z=0 



-R 



dip(R, z) 
dR 
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4 Composite models 



Two families of the models considered in sections 2, 3 give us favorable opportunities for modeling 
the various galaxies and star clusters. But we get more flexible instrument if we use superposition 
of the models of these families. We define the set of parameters as 

p n = {e n , fi n , a n , K n , A\ n) } , n = 1, N. 

here i — number of the selected expression in f ill I) . Then N-component model of the gravitational 
field is the weighted sum of the potentials (J3j), (J4j) 

N N 
n=l n=l 

where w n are assigned normalized weight coefficients. Mass densities are summed up in just the 
same way. For example, if we take N = 2, E\ = 0, > 0, = A^' = 0, and all the other model 
parameters in both components we retain unchanged then we have five free parameters, including 
W\. In such a model force relief field might involve crater either do not involve or involve crater with 
central peak. 



5 Investigation and discussion of orbits 

5.0. We solve Cauchy problem to construct the trajectory of the star. The system of ordinary 
differential equations of the second order is 

r = V^(0 • 

Here r is a position vector of the star with respect to the center. Initial conditions are given for 
arbitrary moment of time to'. 

r(t ) = r , v(£ ) = v , 

where v = r — velocity vector of the star. 

We use dimensionless time (and all the other functions and variables). Interval of integration is 
characterized by number of crossings the plane z = or y = 0. 

System of ordinary differential equations of the second order is transformed to the form of equa- 
tions of the first order in a usual way and is solved numerically. We used Merson's method and 
Dormand-Prince one with variable time step (Hairer, Nersett, Wanner 1990) in order to control 
specified accuracy. We have chosen these methods, as they supplied rather small relative error for 
the integral of energy. 
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5.1. Energy E and ^-component of the angular momentum L z are the integrals of motion for the 
axisymmetrical model (Ogorodnikov 1958) 



E = V (0-^v 2 m + v 2 e ), L z = Rv e , v 2 n 



m - V R 




Here vr, vq, v z are the velocity components in cylindrical coordinates, v m is a meridional velocity, 
E is the total energy with the opposite sign. Using z-component of the angular momentum one can 
write the equations of motion as a system of five equations of the first order (Kutuzov & Ossipkov 



Fr, F z are the components of force F per unit mass. 

We calculate trajectories of the stars either with the initial zero full velocity or with the initial 
zero meridional velocity. That stars we call falling ones. Some periodic orbits were found (Fig. 1). 
Fig. la shows the orbit of the star, falling from the contour of zero-velocity curve. Fig. lb presents 
the orbit of the falling star, with the initial zero full velocity. 

If we assign small perturbations in the initial phase coordinates these orbits conserve the form, 
filling the tube with varying width, which includes initial periodic parent orbit (Fig. 2a). We form 
perturbations such a way so to make star to begin its motion from the equipotential passing through 
the initial point of the parent orbit. As we see in Fig. 2b the orbit does not conserve its form and 
does not include parent orbit if perturbations are rather large. We have similar result while varying 
velocity components. Notice that trajectory does not reach equipotential curve in the half-plane 
z < in Fig. lb, 2a . 

5.2. We consider orbits in triaxial model further. Trajectories of a star, falling from the same 
initial point in the fields with various triaxiality parameters, are shown in Fig. 3. 

It's also interesting to consider trajectory when the inial point lies on the x-axis, while initial 
velocity is directed along the |/-axis. For the axially symmetric model the trajectory lies inside the 
area bounded by circular orbit if the initial velocity (its nonzero component v y ) is smaller than 
circular velocity at the initial point. It lies outside if the initial velocity is larger than circular one. 
The same situation is observed for the triaxial model (Fig. 4) in respect of the quasi-circular velocity 



It's convenient to study orbits in the co-moving meridional plane for the axisymmetrical poten- 
tials (Fig. 5a). We construct orbits in the analogical plane for the triaxial model (Fig. 5b). Bounds 
of the orbit become more "dishevelled" . 

By analogy with Poincare surfaces of section (Binney & Tremaine 1987) we can construct six 
diagrams. The first pair of them has axes x, x. Points with these coordinates are plotted at the 



1981): 




(13) 



(USD- 
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Figure 1. Periodic orbits in axisymmetrical model: the trajectory of the falling star inside the area 
bounded by the zero-velocity curve. 




Figure 2. Effect of varying initial conditions for periodic orbit in Fig. lb. 
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Initial point (marked by dot): 
Xg = 0.5 
Bo = 0.5 

z = 

%o = o 
%> = o 
%a = o 



N„ = 150 



(a) 




Potential KM 
3-axes equipotentials 
Parameters: 
A 2 = 0.01 
a = 2 
, J= 1 
e = 0.5 
i « = 1 



Initial point (marked by dot): 
£0 = 0.5 
2/0 = 0.5 

z = 
v x o = 
ti„o = 
v z o = 



AT S = 150 



Potential KM 
3-axes equipotentials 
Parameters: 
A 2 = l 
a = 2 
»=1 
e = 0.5 

\ K = 1 




Figure 3. Variations of the coefficient A 2 in triaxility parameter r. 




Figure 4. Loop and box orbits in triaxial model 
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Figure 5. Two trajectories in their co-moving planes, (a) — in axysimmetrical model; (b) — in 
triaxial model. Iniatial positions and velocity components are the same for both trajectories 



moments when y — 0, y > or z = 0, z > respectively. Other diagrams could be obtained by 
cyclic permutation of coordinates. Surfaces of section (x, v x ) when y — 0, v y > for the triaxial 
orbits are shown in Fig. 6. Calculations are made for the models with various values of e. We remind 
that the value e — 1 means a spherical model and e = means a model with embedded sharp disk. 
Obviously a stochasticity grows with a flattening. 

We plan to continue the calculation and analysis of orbits in various orbits, describing real galaxies. 
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Figure 6. Surfaces of section for orbits in triaxial models with various flattening. Initial positions 
are the same for all orbits, values of parameters are the same too, except e. N y is the number of 
crossing the plane y = with positive v y . 
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